调节效应(Moderation Effect)是统计学和社会科学中的一个概念,指的是某个第三变量(称为调节变量或调节因子)能够改变两个其他变量之间的关系的强度或方向所产生的作用。通过调节效应的分析,可以更好地理解复杂的变量关系,为研究者提供更加精细化的理论指导和建议。
当一个调节变量存在时,它会通过影响自变量与因变量之间的关系来发挥作用。例如:
调节效应通常用交互图(Interaction Plot)展示。
一项关于气候变化与灾害Chapman and Lickel 2015的研究中,研究者向211名实验参与者讲述非洲发生干旱造成人道主义危机,告诉其中一半的参与者气候变化是造成干旱的原因,另一半参与者未被告知任何关于干旱的原因。接着通过一系列问题让实验参与者评价拒绝援助的正当性,以及了解他们对气候变化的怀疑程度。最后了解参与者捐款的意愿。 实验的目的是了解框架(frame),即是否告知干旱是由于气候变化产生的,对捐助意愿(donate)的影响数据。
假定对气候变化怀疑(skeptic)程度较高的人框架对捐助意愿的效应也比较小,会存在调节效应。
disaster <- read.csv("disaster.csv", header = T)
# 调节效应
moderafit <- lm(donate ~ frame + skeptic + frame:skeptic, data = disaster)
summary(moderafit)
##
## Call:
## lm(formula = donate ~ frame + skeptic + frame:skeptic, data = disaster)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8341 -0.7077 0.1659 0.9101 2.6682
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.02947 0.22632 22.223 <2e-16 ***
## frame 0.67930 0.33091 2.053 0.0413 *
## skeptic -0.13953 0.05790 -2.410 0.0168 *
## frame:skeptic -0.17071 0.08393 -2.034 0.0432 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.234 on 207 degrees of freedom
## Multiple R-squared: 0.1343, Adjusted R-squared: 0.1218
## F-statistic: 10.71 on 3 and 207 DF, p-value: 1.424e-06
中介效应(Mediation Effect)是统计学和社会科学中的一个概念,指一个中介变量(Mediating Variable)在自变量(Independent Variable, \(X\))和因变量(Dependent Variable, \(Y\))之间起中介作用,即部分或全部传递自变量对因变量的影响。中介效应主要用于探讨变量之间的作用机制,揭示因果路径。中介效应的分析能够揭示隐藏的机制和作用路径,为理论研究提供深入的解释,并为实践应用提供精细化的决策依据。
如果 \(X\) 通过一个中介变量 \(M\) 影响 \(Y\),这表明 \(M\) 是 \(X\) 和 \(Y\) 之间的中介变量。公式化表达为: \[ X \rightarrow M \rightarrow Y \] 其中:
中介效应说明 \(X\) 不仅直接影响 \(Y\)(直接效应),还通过 \(M\) 间接影响 \(Y\)(间接效应)。
经典的中介效应模型可以分为三步:
间接效应由 \(a \times b\) 表示,总效应可以分解为: \[ c = c' + (a \times b) \]
上面气候变化与灾害的研究中,框架可能是通过影响正当性,然后正当性再影响捐助意愿的,即中介效应。
library(lavaan)
## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.
# 中介效应
mediamodel <- ' # 直接效应
donate ~ c*frame
# 中介效应
justify ~ a*frame
donate ~ b*justify
# 间接效应 (a*b)
ab := a*b
# 总效应
total := c + (a*b)
'
mediafit <- sem(mediamodel, data = disaster)
summary(mediafit)
## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 211
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## donate ~
## frame (c) 0.212 0.135 1.576 0.115
## justify ~
## frame (a) 0.134 0.127 1.054 0.292
## donate ~
## justify (b) -0.953 0.072 -13.159 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .donate 0.948 0.092 10.271 0.000
## .justify 0.856 0.083 10.271 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## ab -0.128 0.122 -1.051 0.293
## total 0.084 0.181 0.463 0.643
绘制图形
library("semPlot")
semPaths(mediafit,whatLabels = 'est',residuals = F, nCharNodes=0, sizeMan = 12,edge.label.cex = 1.5)
条件过程模型(Conditional Process Model)是一个综合框架,用于同时分析中介效应(Mediation Effect)和调节效应(Moderation Effect)。它探讨自变量通过中介变量影响因变量的机制,同时考虑调节变量如何影响这一机制的强度或方向。条件过程模型是一种强大的分析工具,能够揭示复杂的变量关系,为研究者提供深入的理论解释和实践指导。
条件过程模型结合了中介效应和调节效应,回答以下关键问题:
公式化表示为: \[ Y = b_1 M + b_2 X + b_3 W + b_4 (M \cdot W) + \epsilon \] 其中:
示例:环保政策严格性 (\(X\)) 通过企业创新 (\(M\)) 改善环境质量 (\(Y\)),但监管强度 (\(W\)) 调节这一过程。
示例:政策执行力度 (\(X\)) 和地方经济发展水平 (\(W\)) 的交互作用通过基层治理能力 (\(M\)) 间接影响社会稳定 (\(Y\))。
解释:同时存在调节的中介效应和中介的调节效应,展示了变量之间更复杂的关系。 - \(X → M → Y\):中介效应 - \(W\):调节变量同时作用于 \(X → M\) 和 \(M → Y\) 的路径。
示例:公共政策透明度 (\(X\)) 通过公众信任 (\(M\)) 增强公众满意度 (\(Y\)),但这种过程因社会参与度 (\(W\)) 的不同而改变。
综合调节与中介效应,框架对捐助意愿的直接和间接效应是受到怀疑程度的调节的,即被调节的中介效应。
# 具有调节中介效应的条件过程模型
cpmodel <- ' # 直接效应
donate ~ c1*frame + c2*skeptic + c3*skeptic:frame
# 中介效应
justify ~ a1*frame + a2*skeptic + a3*skeptic:frame
donate ~ b*justify
# 间接效应取决于skeptic的值,需要用平均值(或其他代表值)带入计算间接效应
# 间接效应 (a*b) skeptic取平均值3.38
a1b := (a1+a3*3.38)*b
# 总效应
total := c1 + (a1+a3*3.38)*b
'
cpfit <- sem(cpmodel, data = disaster)
summary(cpfit)
## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 211
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## donate ~
## frame (c1) 0.160 0.264 0.606 0.544
## skeptic (c2) -0.043 0.046 -0.918 0.359
## skptc:frm (c3) 0.015 0.068 0.219 0.827
## justify ~
## frame (a1) -0.562 0.216 -2.606 0.009
## skeptic (a2) 0.105 0.038 2.782 0.005
## skptc:frm (a3) 0.201 0.055 3.675 0.000
## donate ~
## justify (b) -0.923 0.083 -11.113 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .donate 0.943 0.092 10.271 0.000
## .justify 0.648 0.063 10.271 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## a1b -0.108 0.103 -1.054 0.292
## total 0.052 0.285 0.182 0.856
图形
semPaths(cpfit,whatLabels = 'est', layout = "spring",,residuals = F, nCharNodes=0, sizeMan = 8,edge.label.cex = 1)
\[\rho_c=\frac{(\Sigma\lambda_i)^2}{(\Sigma\lambda^2_i + \Sigma\Theta_{ii})}\]
\[\rho_{\nu}=\frac{\Sigma\lambda^2_i}{(\Sigma\lambda^2_i + \Sigma\Theta_{ii})}=\frac{\Sigma\lambda^2_i}{n}\]
案例来自《组织创新气氛量表》(邱皓政,1999),样本是384位企业员工,量表为Likert式6点度量的自陈量表,基于理论和文献先界定影响组织气氛知觉的因素包括组织价值、工作方式、团队合作、领导风格、学习成长、环境气氛等6个因素,每个因素采用3个题项测量,共有18个题项,题项描述性统计如下。
library(haven)
library(tidyverse)
library(modelsummary)
library(lavaan)
library(semPlot)
dat <- read_sav("ch05.sav")
names(dat) <- c("A1","A2","A3","B1","B2","B3","C1","C2","C3","D1","D2","D3","E1","E2","E3","F1","F2","F3")
dat <- zap_labels(dat)
dim(dat)
## [1] 313 18
datasummary_skim(dat, fun_numeric = list('取值'=NUnique, '缺失值'=PercentMissing, '均值'=Mean, '中位数'=Median,'标准差'=SD, '最小值'=Min, '最大值'=Max), output = "data.frame", type = "numeric", fmt_sprintf("%.2f")) |> knitr::kable()
| 取值 | 缺失值 | 均值 | 中位数 | 标准差 | 最小值 | 最大值 | |
|---|---|---|---|---|---|---|---|
| A1 | 6 | 0 | 4.42 | 5.00 | 0.98 | 1.00 | 6.00 |
| A2 | 6 | 0 | 4.31 | 4.00 | 1.02 | 1.00 | 6.00 |
| A3 | 6 | 0 | 4.07 | 4.00 | 0.97 | 1.00 | 6.00 |
| B1 | 6 | 0 | 4.02 | 4.00 | 1.16 | 1.00 | 6.00 |
| B2 | 6 | 0 | 4.25 | 4.00 | 1.16 | 1.00 | 6.00 |
| B3 | 6 | 0 | 4.24 | 4.00 | 1.09 | 1.00 | 6.00 |
| C1 | 6 | 0 | 4.37 | 4.00 | 0.98 | 1.00 | 6.00 |
| C2 | 6 | 0 | 4.34 | 4.00 | 1.03 | 1.00 | 6.00 |
| C3 | 6 | 0 | 4.31 | 4.00 | 1.05 | 1.00 | 6.00 |
| D1 | 6 | 0 | 4.83 | 5.00 | 0.94 | 1.00 | 6.00 |
| D2 | 5 | 0 | 4.95 | 5.00 | 0.84 | 2.00 | 6.00 |
| D3 | 5 | 0 | 4.83 | 5.00 | 0.91 | 2.00 | 6.00 |
| E1 | 6 | 0 | 4.63 | 5.00 | 0.97 | 1.00 | 6.00 |
| E2 | 6 | 0 | 4.73 | 5.00 | 1.01 | 1.00 | 6.00 |
| E3 | 6 | 0 | 4.70 | 5.00 | 0.98 | 1.00 | 6.00 |
| F1 | 6 | 0 | 4.23 | 4.00 | 1.17 | 1.00 | 6.00 |
| F2 | 6 | 0 | 4.63 | 5.00 | 1.09 | 1.00 | 6.00 |
| F3 | 6 | 0 | 4.49 | 5.00 | 0.94 | 1.00 | 6.00 |
lavaan包的模型设定默认不同潜变量之间具有相关性(即允许因子共变),如果需要不同的设定,可以在cfa函数中用\(orthogonal=T\)将协方差约束为0;cfa函数默认是将第1个因子载荷设定为1,如需固定潜变量方差为1,可以通过\(std.lv=TRUE\)来设定。
cfa_model <-'
#定义测量模型
FA =~ L11*A1 + L21*A2 + L31*A3
FB =~ L12*B1 + L22*B2 + L32*B3
FC =~ L13*C1 + L23*C2 + L33*C3
FD =~ L14*D1 + L24*D2 + L34*D3
FE =~ L15*E1 + L25*E2 + L35*E3
FF =~ L16*F1 + L26*F2 + L36*F3'
cfa_fit <- cfa(model = cfa_model,
data = dat,
std.lv = TRUE) # 根据模型设定,固定潜变量方差为1
summary(cfa_fit,
fit.measures = T, # 输出拟合指标
standard = T) # 输出标准化解
## lavaan 0.6-19 ended normally after 31 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 51
##
## Number of observations 313
##
## Model Test User Model:
##
## Test statistic 241.755
## Degrees of freedom 120
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 2842.819
## Degrees of freedom 153
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.955
## Tucker-Lewis Index (TLI) 0.942
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -6751.785
## Loglikelihood unrestricted model (H1) -6630.907
##
## Akaike (AIC) 13605.569
## Bayesian (BIC) 13796.626
## Sample-size adjusted Bayesian (SABIC) 13634.870
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.057
## 90 Percent confidence interval - lower 0.047
## 90 Percent confidence interval - upper 0.067
## P-value H_0: RMSEA <= 0.050 0.132
## P-value H_0: RMSEA >= 0.080 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.052
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## FA =~
## A1 (L11) 0.815 0.051 15.967 0.000 0.815 0.830
## A2 (L21) 0.706 0.055 12.733 0.000 0.706 0.692
## A3 (L31) 0.614 0.054 11.430 0.000 0.614 0.634
## FB =~
## B1 (L12) 0.789 0.062 12.759 0.000 0.789 0.682
## B2 (L22) 0.961 0.058 16.591 0.000 0.961 0.833
## B3 (L32) 0.856 0.056 15.406 0.000 0.856 0.788
## FC =~
## C1 (L13) 0.699 0.053 13.178 0.000 0.699 0.717
## C2 (L23) 0.736 0.056 13.119 0.000 0.736 0.715
## C3 (L33) 0.696 0.058 11.950 0.000 0.696 0.663
## FD =~
## D1 (L14) 0.810 0.045 18.177 0.000 0.810 0.867
## D2 (L24) 0.741 0.040 18.750 0.000 0.741 0.886
## D3 (L34) 0.655 0.047 14.086 0.000 0.655 0.720
## FE =~
## E1 (L15) 0.806 0.046 17.392 0.000 0.806 0.830
## E2 (L25) 0.912 0.046 19.851 0.000 0.912 0.906
## E3 (L35) 0.791 0.047 16.810 0.000 0.791 0.811
## FF =~
## F1 (L16) 0.641 0.066 9.663 0.000 0.641 0.550
## F2 (L26) 0.828 0.058 14.296 0.000 0.828 0.758
## F3 (L36) 0.786 0.049 16.142 0.000 0.786 0.837
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## FA ~~
## FB 0.542 0.054 9.989 0.000 0.542 0.542
## FC 0.494 0.061 8.102 0.000 0.494 0.494
## FD 0.417 0.058 7.176 0.000 0.417 0.417
## FE 0.526 0.052 10.080 0.000 0.526 0.526
## FF 0.695 0.046 15.148 0.000 0.695 0.695
## FB ~~
## FC 0.697 0.047 14.906 0.000 0.697 0.697
## FD 0.447 0.055 8.127 0.000 0.447 0.447
## FE 0.575 0.048 12.091 0.000 0.575 0.575
## FF 0.391 0.061 6.403 0.000 0.391 0.391
## FC ~~
## FD 0.522 0.055 9.493 0.000 0.522 0.522
## FE 0.603 0.050 12.127 0.000 0.603 0.603
## FF 0.600 0.054 11.032 0.000 0.600 0.600
## FD ~~
## FE 0.557 0.046 12.046 0.000 0.557 0.557
## FF 0.316 0.062 5.137 0.000 0.316 0.316
## FE ~~
## FF 0.443 0.056 7.944 0.000 0.443 0.443
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .A1 0.300 0.046 6.509 0.000 0.300 0.311
## .A2 0.544 0.055 9.980 0.000 0.544 0.522
## .A3 0.561 0.052 10.700 0.000 0.561 0.598
## .B1 0.716 0.068 10.525 0.000 0.716 0.535
## .B2 0.408 0.057 7.211 0.000 0.408 0.306
## .B3 0.447 0.052 8.602 0.000 0.447 0.379
## .C1 0.461 0.050 9.305 0.000 0.461 0.485
## .C2 0.519 0.055 9.354 0.000 0.519 0.489
## .C3 0.618 0.061 10.173 0.000 0.618 0.560
## .D1 0.217 0.031 7.055 0.000 0.217 0.248
## .D2 0.151 0.024 6.204 0.000 0.151 0.215
## .D3 0.399 0.037 10.799 0.000 0.399 0.482
## .E1 0.292 0.032 9.243 0.000 0.292 0.311
## .E2 0.181 0.030 6.113 0.000 0.181 0.179
## .E3 0.325 0.033 9.733 0.000 0.325 0.342
## .F1 0.948 0.083 11.400 0.000 0.948 0.698
## .F2 0.507 0.058 8.700 0.000 0.507 0.425
## .F3 0.265 0.042 6.276 0.000 0.265 0.300
## FA 1.000 1.000 1.000
## FB 1.000 1.000 1.000
## FC 1.000 1.000 1.000
## FD 1.000 1.000 1.000
## FE 1.000 1.000 1.000
## FF 1.000 1.000 1.000
summary报告第一部分为模型基本信息,估计方法为最大似然估计,优化方法为非线性最小优化,待估参数51个,样本量313。
summary报告第二部分为拟合指标,包括卡方值和检验结果、CFI、RMSEA、SRMR等。
fitmeasures(cfa_fit, c("chisq", "df", "pvalue", "nfi","nnfi","cfi", "rmsea", "srmr"))
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 241.755 120.000 0.000 0.915 0.942 0.955 0.057 0.052
summary报告第三部分为参数估计结果,包括因子载荷、潜变量协方差、测量误差和潜变量的方差。
std.lv对应的是潜变量被标准化时的载荷,std.all对应的是所有变量都被标准化时的载荷(类似于标准化的回归系数)。模型拟合后可以依据拟合协方差矩阵列出估计的测量变量的方差和协方差,这些模型导出数与实际观察值之间的差距即为残差。通过分析标准化残差的大小和分布,可以了解测量变量间关系的拟合程度,确定存在问题的题项。
fitted(cfa_fit) # 提取拟合协方差矩阵
## $cov
## A1 A2 A3 B1 B2 B3 C1 C2 C3 D1 D2 D3
## A1 0.965
## A2 0.576 1.043
## A3 0.500 0.433 0.937
## B1 0.349 0.302 0.262 1.338
## B2 0.425 0.368 0.320 0.759 1.332
## B3 0.378 0.328 0.285 0.675 0.823 1.179
## C1 0.282 0.244 0.212 0.385 0.469 0.417 0.950
## C2 0.297 0.257 0.223 0.405 0.494 0.439 0.515 1.061
## C3 0.280 0.243 0.211 0.383 0.467 0.415 0.487 0.512 1.102
## D1 0.275 0.239 0.207 0.286 0.348 0.310 0.296 0.312 0.295 0.873
## D2 0.252 0.218 0.190 0.262 0.319 0.284 0.271 0.285 0.269 0.601 0.700
## D3 0.223 0.193 0.168 0.231 0.282 0.251 0.239 0.252 0.238 0.531 0.486 0.829
## E1 0.345 0.299 0.260 0.365 0.445 0.396 0.340 0.358 0.338 0.363 0.332 0.294
## E2 0.391 0.339 0.294 0.413 0.504 0.448 0.385 0.405 0.383 0.411 0.376 0.333
## E3 0.339 0.294 0.255 0.359 0.437 0.389 0.334 0.352 0.332 0.357 0.326 0.289
## F1 0.363 0.315 0.274 0.198 0.241 0.214 0.269 0.283 0.268 0.164 0.150 0.133
## F2 0.469 0.407 0.353 0.255 0.311 0.277 0.347 0.366 0.346 0.212 0.194 0.171
## F3 0.446 0.386 0.335 0.243 0.296 0.263 0.330 0.347 0.328 0.201 0.184 0.163
## E1 E2 E3 F1 F2 F3
## A1
## A2
## A3
## B1
## B2
## B3
## C1
## C2
## C3
## D1
## D2
## D3
## E1 0.942
## E2 0.735 1.013
## E3 0.638 0.721 0.951
## F1 0.229 0.259 0.225 1.359
## F2 0.296 0.335 0.291 0.531 1.193
## F3 0.281 0.318 0.276 0.504 0.651 0.883
cfares <- resid(cfa_fit) # 提取残差
cfares$cov
## A1 A2 A3 B1 B2 B3 C1 C2 C3 D1 D2
## A1 0.000
## A2 -0.003 0.000
## A3 0.022 -0.037 0.000
## B1 -0.080 0.057 0.037 0.000
## B2 -0.009 0.104 -0.057 0.001 0.000
## B3 -0.021 0.053 -0.011 0.002 -0.001 0.000
## C1 -0.085 0.022 -0.054 0.022 0.015 0.024 0.000
## C2 -0.019 0.039 -0.101 0.005 -0.020 -0.041 0.042 0.000
## C3 0.070 0.167 0.079 -0.004 0.002 0.006 -0.050 -0.004 0.000
## D1 0.030 0.113 -0.048 -0.012 0.026 0.043 -0.052 -0.006 0.022 0.000
## D2 -0.033 0.021 -0.071 -0.031 -0.041 -0.029 -0.076 -0.003 0.018 0.008 0.000
## D3 -0.010 0.075 -0.050 -0.005 0.057 0.067 0.094 0.114 0.117 -0.028 0.004
## E1 -0.017 0.015 0.087 0.021 0.073 0.063 -0.001 -0.077 0.047 0.055 0.013
## E2 -0.047 0.100 0.032 -0.029 -0.034 -0.024 0.024 -0.034 0.016 -0.021 -0.054
## E3 -0.103 0.105 -0.010 -0.031 -0.013 0.019 0.005 -0.004 0.045 0.005 -0.019
## F1 -0.018 -0.060 0.058 -0.035 0.097 0.021 0.070 0.118 0.111 0.107 0.130
## F2 -0.007 -0.043 -0.064 -0.122 -0.060 0.004 -0.080 -0.048 0.121 -0.015 -0.005
## F3 0.031 -0.007 0.018 -0.027 0.040 0.015 -0.006 -0.057 0.024 -0.033 -0.044
## D3 E1 E2 E3 F1 F2 F3
## A1
## A2
## A3
## B1
## B2
## B3
## C1
## C2
## C3
## D1
## D2
## D3 0.000
## E1 0.083 0.000
## E2 0.074 -0.001 0.000
## E3 0.045 -0.016 0.010 0.000
## F1 0.187 0.076 0.109 0.097 0.000
## F2 0.038 -0.056 -0.012 -0.019 0.019 0.000
## F3 0.028 -0.017 0.001 -0.016 -0.032 0.009 0.000
range(cfares$cov); median(cfares$cov) # 残差的全距和中位数
## [1] -0.1216818 0.1866437
## [1] 7.123659e-07
resid(cfa_fit, type="standardized") # 提取标准化残差,方便与标准正态分布比较
## $type
## [1] "standardized"
##
## $cov
## A1 A2 A3 B1 B2 B3 C1 C2 C3 D1 D2
## A1 0.000
## A2 -0.363 0.000
## A3 1.765 -1.670 0.000
## B1 -2.017 1.218 0.774 0.000
## B2 -0.326 2.626 -1.426 0.055 0.000
## B3 -0.700 1.340 -0.288 0.081 -0.146 0.000
## C1 -2.886 0.559 -1.395 0.552 0.522 0.784 0.000
## C2 -0.611 0.950 -2.455 0.117 -0.662 -1.305 2.421 0.000
## C3 1.876 3.832 1.833 -0.099 0.044 0.162 -2.676 -0.185 0.000
## D1 1.262 3.253 -1.394 -0.301 0.901 1.417 -1.900 -0.216 0.653 0.000
## D2 -1.688 0.684 -2.317 -0.888 -1.700 -1.090 -3.170 -0.104 0.599 3.454 0.000
## D3 -0.301 1.933 -1.283 -0.112 1.480 1.767 2.707 3.162 3.024 -4.148 0.748
## E1 -0.643 0.417 2.439 0.528 2.335 1.975 -0.045 -2.494 1.359 2.226 0.609
## E2 -2.361 3.078 0.956 -0.781 -1.407 -0.873 0.943 -1.297 0.479 -1.071 -3.481
## E3 -3.617 2.895 -0.260 -0.760 -0.408 0.568 0.172 -0.119 1.255 0.213 -0.859
## F1 -0.456 -1.244 1.205 -0.544 1.656 0.373 1.449 2.317 2.177 2.172 2.947
## F2 -0.295 -1.216 -1.743 -2.351 -1.483 0.104 -2.284 -1.330 2.968 -0.463 -0.189
## F3 1.823 -0.259 0.661 -0.666 1.394 0.490 -0.222 -2.153 0.745 -1.412 -2.229
## D3 E1 E2 E3 F1 F2 F3
## A1
## A2
## A3
## B1
## B2
## B3
## C1
## C2
## C3
## D1
## D2
## D3 0.000
## E1 2.629 0.000
## E2 2.415 -0.342 0.000
## E3 1.378 -1.831 2.050 0.000
## F1 3.664 1.520 2.231 1.915 0.000
## F2 0.933 -1.571 -0.403 -0.528 0.649 0.000
## F3 0.835 -0.641 0.056 -0.579 -2.117 1.473 0.000
cfastd <- cfares$cov[lower.tri(unclass(cfares$cov),diag = F)]
ggplot(mapping = aes(sample=cfastd)) +
geom_qq() +
geom_qq_line() # 绘制QQ残差散点图分析是否服从正态分布
modificationindices(cfa_fit, sort. = T, maximum.number = 6)
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 96 FC =~ D3 19.838 0.241 0.241 0.265 0.265
## 63 FA =~ C3 14.826 0.289 0.289 0.276 0.276
## 118 FE =~ A1 14.048 -0.258 -0.258 -0.262 -0.262
## 265 D1 ~~ D2 13.890 0.192 0.192 1.061 1.061
## 266 D1 ~~ D3 13.426 -0.132 -0.132 -0.448 -0.448
## 175 A2 ~~ E1 12.247 -0.097 -0.097 -0.243 -0.243
根据内部拟合检验部分提出的各项内部拟合指标,对测量模型进行检验。
cfa_lambda <- inspect(cfa_fit,what="std")$lambda
cfa_lambda <- cfa_lambda[cfa_lambda!=0]
print(cfa_lambda, digits = 2) # 检验因子载荷
## [1] 0.83 0.69 0.63 0.68 0.83 0.79 0.72 0.71 0.66 0.87 0.89 0.72 0.83 0.91 0.81
## [16] 0.55 0.76 0.84
library(semTools)
compRelSEM(cfa_fit) # 计算组合信度CR
## FA FB FC FD FE FF
## 0.769 0.812 0.743 0.869 0.889 0.748
AVE(cfa_fit) # 计算平均变异萃取量AVE
## FA FB FC FD FE FF
## 0.523 0.592 0.487 0.681 0.725 0.499
summary报告第三部分的潜变量协方差,即相关系数)。cfa_res <- summary(cfa_fit, standard=T)
cfa_res$pe[43:57,] |>
select(lhs, op, rhs, est, se, z, std.all, pvalue) |>
mutate(across(where(is.numeric), ~round(., 3)))
## lhs op rhs est se z std.all pvalue
## 43 FA ~~ FB 0.542 0.054 9.989 0.542 0
## 44 FA ~~ FC 0.494 0.061 8.102 0.494 0
## 45 FA ~~ FD 0.417 0.058 7.176 0.417 0
## 46 FA ~~ FE 0.526 0.052 10.080 0.526 0
## 47 FA ~~ FF 0.695 0.046 15.148 0.695 0
## 48 FB ~~ FC 0.697 0.047 14.906 0.697 0
## 49 FB ~~ FD 0.447 0.055 8.127 0.447 0
## 50 FB ~~ FE 0.575 0.048 12.091 0.575 0
## 51 FB ~~ FF 0.391 0.061 6.403 0.391 0
## 52 FC ~~ FD 0.522 0.055 9.493 0.522 0
## 53 FC ~~ FE 0.603 0.050 12.127 0.603 0
## 54 FC ~~ FF 0.600 0.054 11.032 0.600 0
## 55 FD ~~ FE 0.557 0.046 12.046 0.557 0
## 56 FD ~~ FF 0.316 0.062 5.137 0.316 0
## 57 FE ~~ FF 0.443 0.056 7.944 0.443 0
根据前面CFA分析检验的结果,MI指数显示测量变量C3与潜变量A之间的关系如果纳入模型能够提升拟合度14.82,增加因子载荷0.28,综合考量是贡献最大的参数改善。因此需要调整假设和设定,测量变量C3同时受到潜变量C和A的影响,其他参数设定不变。
new_model <-'
#define the measurement model
FA =~ L11*A1 + L21*A2 + L31*A3 + L37 * C3
FB =~ L12*B1 + L22*B2 + L32*B3
FC =~ L13*C1 + L23*C2 + L33*C3
FD =~ L14*D1 + L24*D2 + L34*D3
FE =~ L15*E1 + L25*E2 + L35*E3
FF =~ L16*F1 + L26*F2 + L36*F3'
new_fit <- cfa(new_model,
data = dat,
std.lv = TRUE)
summary(new_fit)
## lavaan 0.6-19 ended normally after 30 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 52
##
## Number of observations 313
##
## Model Test User Model:
##
## Test statistic 226.794
## Degrees of freedom 119
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## FA =~
## A1 (L11) 0.805 0.051 15.803 0.000
## A2 (L21) 0.711 0.055 12.859 0.000
## A3 (L31) 0.616 0.054 11.510 0.000
## C3 (L37) 0.268 0.065 4.128 0.000
## FB =~
## B1 (L12) 0.789 0.062 12.766 0.000
## B2 (L22) 0.961 0.058 16.591 0.000
## B3 (L32) 0.855 0.056 15.407 0.000
## FC =~
## C1 (L13) 0.729 0.054 13.605 0.000
## C2 (L23) 0.755 0.057 13.284 0.000
## C3 (L33) 0.533 0.066 8.080 0.000
## FD =~
## D1 (L14) 0.811 0.045 18.185 0.000
## D2 (L24) 0.741 0.040 18.743 0.000
## D3 (L34) 0.655 0.047 14.078 0.000
## FE =~
## E1 (L15) 0.806 0.046 17.389 0.000
## E2 (L25) 0.912 0.046 19.857 0.000
## E3 (L35) 0.791 0.047 16.809 0.000
## FF =~
## F1 (L16) 0.641 0.066 9.656 0.000
## F2 (L26) 0.828 0.058 14.305 0.000
## F3 (L36) 0.786 0.049 16.137 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## FA ~~
## FB 0.545 0.054 10.046 0.000
## FC 0.409 0.067 6.078 0.000
## FD 0.424 0.058 7.319 0.000
## FE 0.533 0.052 10.283 0.000
## FF 0.700 0.046 15.368 0.000
## FB ~~
## FC 0.668 0.049 13.556 0.000
## FD 0.447 0.055 8.128 0.000
## FE 0.574 0.048 12.089 0.000
## FF 0.391 0.061 6.401 0.000
## FC ~~
## FD 0.491 0.057 8.614 0.000
## FE 0.569 0.052 10.857 0.000
## FF 0.545 0.059 9.309 0.000
## FD ~~
## FE 0.557 0.046 12.045 0.000
## FF 0.316 0.062 5.137 0.000
## FE ~~
## FF 0.443 0.056 7.943 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .A1 0.317 0.045 6.981 0.000
## .A2 0.538 0.054 9.962 0.000
## .A3 0.557 0.052 10.694 0.000
## .C3 0.629 0.057 10.951 0.000
## .B1 0.715 0.068 10.524 0.000
## .B2 0.408 0.057 7.221 0.000
## .B3 0.447 0.052 8.608 0.000
## .C1 0.418 0.051 8.222 0.000
## .C2 0.491 0.057 8.612 0.000
## .D1 0.216 0.031 7.036 0.000
## .D2 0.151 0.024 6.209 0.000
## .D3 0.400 0.037 10.802 0.000
## .E1 0.293 0.032 9.249 0.000
## .E2 0.181 0.030 6.108 0.000
## .E3 0.325 0.033 9.736 0.000
## .F1 0.949 0.083 11.402 0.000
## .F2 0.506 0.058 8.689 0.000
## .F3 0.265 0.042 6.283 0.000
## FA 1.000
## FB 1.000
## FC 1.000
## FD 1.000
## FE 1.000
## FF 1.000
fitmeasures(cfa_fit, c("chisq", "df", "pvalue", "nfi","nnfi","cfi", "rmsea", "srmr")) # 原设定模型
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 241.755 120.000 0.000 0.915 0.942 0.955 0.057 0.052
fitmeasures(new_fit, c("chisq", "df", "pvalue", "nfi","nnfi","cfi", "rmsea", "srmr")) # 调整后模型
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 226.794 119.000 0.000 0.920 0.948 0.960 0.054 0.048
semPaths(cfa_fit, what = "std")
semPaths(new_fit, what = "std")
路径分析所包括的变量都是外显观察变量,不涉及潜变量,所用的操作可采用结构方程分析技术,相当于没有测量模型的结构方程。当模型结构过于复杂时,会出现自由度\(df\le0\)的情况,无法计算拟合指数。
研究假设组织气氛知觉包括6个变量:组织价值(VALUE)、工作方式(JOBSTYLE)、团队合作(TEAMWORK)、领导风格(LEADERSH)、学习成长(LEARNING)、环境气氛(ENVIRON)。这6个变量均会影响员工的组织承诺感(commit),进而影响员工的工作绩效(out),此时,组织承诺作为中介变量,同时年资(TENURE)是控制变量,既影响组织承诺,也影响工作绩效。 整个模型有9个观察变量,因此测量数据\(DP=(9\times10)/2=45\)。模型设定6个组织气氛知觉变量之间都相关,因此有15个相关系数参数,7个外生变量有7个方差参数,2个内生变量有2个解释残差(即回归方程的误差项),9个路径回归系数,所以待估参数\(t=15+7+2+9=33\)个。\(t<DP\)因此模型时可识别的。
dat <- read_sav("ch07.sav")
datasummary_skim(dat, fun_numeric = list('取值'=NUnique, '缺失值'=PercentMissing, '均值'=Mean, '中位数'=Median,'标准差'=SD, '最小值'=Min, '最大值'=Max), output = "data.frame", type = "numeric", fmt_sprintf("%.2f")) |> knitr::kable()
| 取值 | 缺失值 | 均值 | 中位数 | 标准差 | 最小值 | 最大值 | |
|---|---|---|---|---|---|---|---|
| out | 9 | 0 | 4.22 | 4.33 | 0.63 | 2.00 | 5.00 |
| commit | 22 | 0 | 9.49 | 9.67 | 1.59 | 4.00 | 12.00 |
| VALUE | 15 | 0 | 4.27 | 4.33 | 0.81 | 1.00 | 6.00 |
| JOBSTYLE | 15 | 0 | 4.18 | 4.33 | 0.95 | 1.00 | 6.00 |
| TEAMWORK | 14 | 0 | 4.33 | 4.33 | 0.83 | 1.33 | 6.00 |
| LEADERSH | 12 | 0 | 4.87 | 5.00 | 0.77 | 2.00 | 6.00 |
| LEARNING | 14 | 0 | 4.71 | 5.00 | 0.88 | 1.67 | 6.00 |
| ENVIRONM | 14 | 0 | 4.43 | 4.33 | 0.88 | 1.00 | 6.00 |
| TENURE | 106 | 0 | 9.12 | 4.00 | 9.23 | 0.08 | 32.00 |
path1 <-'
#part1:设定结构模型
commit~ a1*VALUE+a2*JOBSTYLE+a3*TEAMWORK+a4*LEADERSH
+a5*LEARNING+a6*ENVIRONM+a7*TENURE
out~c*TENURE+b*commit
#part2:设定外生变量的相关系数
VALUE ~~ JOBSTYLE
VALUE ~~ TEAMWORK
VALUE ~~ LEADERSH
VALUE ~~ LEARNING
VALUE ~~ ENVIRONM
JOBSTYLE ~~ TEAMWORK
JOBSTYLE ~~ LEADERSH
JOBSTYLE ~~ LEARNING
JOBSTYLE ~~ ENVIRONM
TEAMWORK ~~ LEADERSH
TEAMWORK ~~ LEARNING
TEAMWORK ~~ ENVIRONM
LEADERSH ~~ LEARNING
LEADERSH ~~ ENVIRONM
LEARNING ~~ ENVIRONM
# 设定相关系数为0
TENURE ~~ 0*VALUE
TENURE ~~ 0*JOBSTYLE
TENURE ~~ 0*TEAMWORK
TENURE ~~ 0*LEADERSH
TENURE ~~ 0*LEARNING
TENURE ~~ 0*ENVIRONM
#part3:定义间接效应 (a*b)
a1b := a1*b
a2b := a2*b
a3b := a3*b
a4b := a4*b
a5b := a5*b
a6b := a6*b
a7b := a7*b
#part4:定义年资TENURE的总效应
total := c + a7b'
summary报告第一部分为模型基本信息,估计方法为最大似然估计,优化方法为非线性最小优化,待估参数33个,样本量281.
summary报告第二部分为拟合度检验,包括卡方值和检验结果、CFI、RMSEA、SRMR等。
- 模型拟合度分析:对照前面的标准,CFI值0.98大于标准0.9,SRMR值0.06小于标准0.08,表明拟合比较理想;RMSEA为0.067略微大于标准0.05,表明模型拟合不理想;综合来看,理论模型拟合度较好。
summary报告第三部分为参数估计结果,包括路径回归系数、外生变量协方差、外生变量方差和内生变量的解释残差、直接效应和间接效应参数估计。
- 回归模型(Regressions)部分为路径回归系数估计值相关结果,std.all对应的是标准化的回归系数。团队合作(TEAMWORK)和领导风格(LEADERSH)对组织承诺(commit)的回归系数不显著。其余回归系数均显著。
- 协方差(Covariances)部分为外生变量协方差,由于模型设定年资(TENURE)与其他外生变量不相关,因此固定为0,std.all是标准化后的协方差,即外生变量的相关系数。协方差估计结果均显著,说明外生变量间相关性存在。
- 方差(Variances)部分为内生变量的解释残差(.commit和.out)和外生变量的方差。
- 定义参数估计(Defined Parameters)部分为总效应和间接效应系数相关信息。同样,团队合作(TEAMWORK)和领导风格(LEADERSH)的间接效应系数a3b和a4b不显著,显然这与回归模型的结论一致。
path1_fit <- sem(path1, data = dat)
path1_res <- summary(path1_fit, fit.measures=T,standard = T)
path1_res
## lavaan 0.6-19 ended normally after 43 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 33
##
## Number of observations 281
##
## Model Test User Model:
##
## Test statistic 27.106
## Degrees of freedom 12
## P-value (Chi-square) 0.007
##
## Model Test Baseline Model:
##
## Test statistic 866.406
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.982
## Tucker-Lewis Index (TLI) 0.945
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3522.025
## Loglikelihood unrestricted model (H1) -3508.473
##
## Akaike (AIC) 7110.051
## Bayesian (BIC) 7230.116
## Sample-size adjusted Bayesian (SABIC) 7125.475
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.067
## 90 Percent confidence interval - lower 0.033
## 90 Percent confidence interval - upper 0.101
## P-value H_0: RMSEA <= 0.050 0.181
## P-value H_0: RMSEA >= 0.080 0.287
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.060
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## commit ~
## VALUE (a1) 0.428 0.116 3.687 0.000 0.428 0.218
## JOBSTYLE (a2) 0.242 0.100 2.409 0.016 0.242 0.146
## TEAMWORK (a3) 0.106 0.118 0.895 0.371 0.106 0.056
## LEADERSH (a4) 0.108 0.119 0.912 0.362 0.108 0.053
## LEARNING (a5) 0.363 0.115 3.159 0.002 0.363 0.203
## ENVIRONM (a6) 0.276 0.106 2.611 0.009 0.276 0.153
## TENURE (a7) 0.022 0.008 2.816 0.005 0.022 0.129
## out ~
## TENURE (c) 0.011 0.004 2.964 0.003 0.011 0.157
## commit (b) 0.172 0.021 8.135 0.000 0.172 0.430
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## VALUE ~~
## JOBSTYLE 0.367 0.051 7.244 0.000 0.367 0.479
## TEAMWORK 0.271 0.043 6.288 0.000 0.271 0.405
## LEADERSH 0.206 0.039 5.274 0.000 0.206 0.331
## LEARNING 0.344 0.047 7.328 0.000 0.344 0.486
## ENVIRONM 0.369 0.047 7.782 0.000 0.369 0.524
## JOBSTYLE ~~
## TEAMWORK 0.420 0.053 7.853 0.000 0.420 0.530
## LEADERSH 0.328 0.048 6.829 0.000 0.328 0.446
## LEARNING 0.461 0.057 8.078 0.000 0.461 0.550
## ENVIRONM 0.269 0.052 5.149 0.000 0.269 0.323
## TEAMWORK ~~
## LEADERSH 0.318 0.043 7.443 0.000 0.318 0.496
## LEARNING 0.390 0.049 7.884 0.000 0.390 0.533
## ENVIRONM 0.368 0.049 7.571 0.000 0.368 0.506
## LEADERSH ~~
## LEARNING 0.386 0.047 8.259 0.000 0.386 0.566
## ENVIRONM 0.249 0.043 5.800 0.000 0.249 0.369
## LEARNING ~~
## ENVIRONM 0.332 0.050 6.626 0.000 0.332 0.430
## VALUE ~~
## TENURE 0.000 0.000 0.000
## JOBSTYLE ~~
## TENURE 0.000 0.000 0.000
## TEAMWORK ~~
## TENURE 0.000 0.000 0.000
## LEADERSH ~~
## TENURE 0.000 0.000 0.000
## LEARNING ~~
## TENURE 0.000 0.000 0.000
## ENVIRONM ~~
## TENURE 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .commit 1.456 0.123 11.853 0.000 1.456 0.586
## .out 0.305 0.026 11.853 0.000 0.305 0.773
## VALUE 0.647 0.055 11.853 0.000 0.647 1.000
## JOBSTYLE 0.907 0.077 11.853 0.000 0.907 1.000
## TEAMWORK 0.691 0.058 11.853 0.000 0.691 1.000
## LEADERSH 0.598 0.050 11.853 0.000 0.598 1.000
## LEARNING 0.776 0.065 11.853 0.000 0.776 1.000
## ENVIRONM 0.765 0.065 11.853 0.000 0.765 1.000
## TENURE 84.834 7.157 11.853 0.000 84.834 1.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## a1b 0.073 0.022 3.358 0.001 0.073 0.094
## a2b 0.041 0.018 2.310 0.021 0.041 0.063
## a3b 0.018 0.020 0.889 0.374 0.018 0.024
## a4b 0.019 0.021 0.906 0.365 0.019 0.023
## a5b 0.062 0.021 2.945 0.003 0.062 0.087
## a6b 0.047 0.019 2.486 0.013 0.047 0.066
## a7b 0.004 0.001 2.661 0.008 0.004 0.055
## total 0.014 0.004 3.787 0.000 0.014 0.212
路径分析的目的之一是探讨内生变量被外生变量解释的总体效应、直接效应和间接效应。尽管采用分阶段回归的方式也能够估计出各个效应值,但是结构方程技术更方便对这些效应进行统计检验。
直接效应:summary报告参数估计的回归模型部分的路径系数便是直接效应。例如组织价值(VALUE)对组织承诺(commit)的直接效应为a1,即0.428,而年资(TENURE)对工作绩效(out)的直接效应为c,即0.011.
间接效应:外生变量通过中介变量对另一个内生变量产生的影响,便是间接效应。例如间接效应系数a1b代表组织价值对组织承诺的效应(a1=0.428)通过组织承诺对工作绩效的直接效应(b=0.172)间接影响到工作绩效,形成间接效应(a1b),并且数值上间接效应等于两段直接效应的乘积,即\(a1b=a1\times b=0.428\times 0.172=0.073\),统计检验显示其p值为0.001,该间接效应具有统计上的意义。
总效应:一个变量对结果变量直接效应和间接效应的总和便是总效应。例如,年资(TENURE)对工作绩效的直接效应(c=0.011),年资通过组织承诺对工作绩效的间接效应(a7b=0.004),那么年资对工作绩效的总效应\(total=c+a7b=0.011+0.004\simeq 0.015\),统计检验其p值为0,该总效应具有统计上的意义。
path1_res$pe[c(1:9,40:47),] |>
select(lhs, op, rhs, est, se, z, std.all, pvalue) |>
mutate(across(where(is.numeric), ~round(., 3)))
## lhs op rhs est se z std.all pvalue
## 1 commit ~ VALUE 0.428 0.116 3.687 0.218 0.000
## 2 commit ~ JOBSTYLE 0.242 0.100 2.409 0.146 0.016
## 3 commit ~ TEAMWORK 0.106 0.118 0.895 0.056 0.371
## 4 commit ~ LEADERSH 0.108 0.119 0.912 0.053 0.362
## 5 commit ~ LEARNING 0.363 0.115 3.159 0.203 0.002
## 6 commit ~ ENVIRONM 0.276 0.106 2.611 0.153 0.009
## 7 commit ~ TENURE 0.022 0.008 2.816 0.129 0.005
## 8 out ~ TENURE 0.011 0.004 2.964 0.157 0.003
## 9 out ~ commit 0.172 0.021 8.135 0.430 0.000
## 40 a1b := a1*b 0.073 0.022 3.358 0.094 0.001
## 41 a2b := a2*b 0.041 0.018 2.310 0.063 0.021
## 42 a3b := a3*b 0.018 0.020 0.889 0.024 0.374
## 43 a4b := a4*b 0.019 0.021 0.906 0.023 0.365
## 44 a5b := a5*b 0.062 0.021 2.945 0.087 0.003
## 45 a6b := a6*b 0.047 0.019 2.486 0.066 0.013
## 46 a7b := a7*b 0.004 0.001 2.661 0.055 0.008
## 47 total := c+a7b 0.014 0.004 3.787 0.212 0.000
modificationindices(path1_fit, sort. = T, maximum.number = 8)
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 48 commit ~~ out 15.228 -0.244 -0.244 -0.366 -0.366
## 63 commit ~ out 15.228 -0.800 -0.800 -0.319 -0.319
## 68 out ~ LEARNING 13.063 0.158 0.158 0.222 0.222
## 67 out ~ LEADERSH 9.170 0.140 0.140 0.173 0.173
## 64 out ~ VALUE 7.031 0.127 0.127 0.162 0.162
## 66 out ~ TEAMWORK 6.771 0.115 0.115 0.152 0.152
## 125 TENURE ~ ENVIRONM 5.352 1.453 1.453 0.138 0.138
## 120 TENURE ~ VALUE 4.581 1.462 1.462 0.128 0.128
先增加一条学习成长到工作绩效的直接效应,设定新的模型,试试看。
path2 <- '#part1:set the structure model
commit~ a1*VALUE+a2*JOBSTYLE+a3*TEAMWORK+a4*LEADERSH
+a5*LEARNING+a6*ENVIRONM+a7*TENURE
out~c*TENURE+b*commit + d*LEARNING # 添加直接效应d
#part2:set the correlation
VALUE ~~ JOBSTYLE
VALUE ~~ TEAMWORK
VALUE ~~ LEADERSH
VALUE ~~ LEARNING
VALUE ~~ ENVIRONM
JOBSTYLE ~~ TEAMWORK
JOBSTYLE ~~ LEADERSH
JOBSTYLE ~~ LEARNING
JOBSTYLE ~~ ENVIRONM
TEAMWORK ~~ LEADERSH
TEAMWORK ~~ LEARNING
TEAMWORK ~~ ENVIRONM
LEADERSH ~~ LEARNING
LEADERSH ~~ ENVIRONM
LEARNING ~~ ENVIRONM
# fix
TENURE ~~ 0*VALUE
TENURE ~~ 0*JOBSTYLE
TENURE ~~ 0*TEAMWORK
TENURE ~~ 0*LEADERSH
TENURE ~~ 0*LEARNING
TENURE ~~ 0*ENVIRONM
#part3:define indirect effect (a*b)
a1b := a1*b
a2b := a2*b
a3b := a3*b
a4b := a4*b
a5b := a5*b
a6b := a6*b
a7b := a7*b
#part4:define total effect of COMMIT
total := c + a7b'
path2_fit <- sem(path2, data = dat)
summary(path2_fit, fit.measures=T, standard=T)
## lavaan 0.6-19 ended normally after 44 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 34
##
## Number of observations 281
##
## Model Test User Model:
##
## Test statistic 13.714
## Degrees of freedom 11
## P-value (Chi-square) 0.249
##
## Model Test Baseline Model:
##
## Test statistic 866.406
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.997
## Tucker-Lewis Index (TLI) 0.989
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3515.330
## Loglikelihood unrestricted model (H1) -3508.473
##
## Akaike (AIC) 7098.659
## Bayesian (BIC) 7222.363
## Sample-size adjusted Bayesian (SABIC) 7114.551
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.030
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.073
## P-value H_0: RMSEA <= 0.050 0.736
## P-value H_0: RMSEA >= 0.080 0.024
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.041
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## commit ~
## VALUE (a1) 0.428 0.116 3.687 0.000 0.428 0.218
## JOBSTYLE (a2) 0.242 0.100 2.409 0.016 0.242 0.146
## TEAMWORK (a3) 0.106 0.118 0.895 0.371 0.106 0.056
## LEADERSH (a4) 0.108 0.119 0.912 0.362 0.108 0.053
## LEARNING (a5) 0.363 0.115 3.159 0.002 0.363 0.203
## ENVIRONM (a6) 0.276 0.106 2.611 0.009 0.276 0.153
## TENURE (a7) 0.022 0.008 2.816 0.005 0.022 0.129
## out ~
## TENURE (c) 0.011 0.004 3.192 0.001 0.011 0.166
## commit (b) 0.125 0.024 5.201 0.000 0.125 0.315
## LEARNING (d) 0.158 0.043 3.706 0.000 0.158 0.222
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## VALUE ~~
## JOBSTYLE 0.367 0.051 7.244 0.000 0.367 0.479
## TEAMWORK 0.271 0.043 6.288 0.000 0.271 0.405
## LEADERSH 0.206 0.039 5.274 0.000 0.206 0.331
## LEARNING 0.344 0.047 7.328 0.000 0.344 0.486
## ENVIRONM 0.369 0.047 7.782 0.000 0.369 0.524
## JOBSTYLE ~~
## TEAMWORK 0.420 0.053 7.853 0.000 0.420 0.530
## LEADERSH 0.328 0.048 6.829 0.000 0.328 0.446
## LEARNING 0.461 0.057 8.078 0.000 0.461 0.550
## ENVIRONM 0.269 0.052 5.149 0.000 0.269 0.323
## TEAMWORK ~~
## LEADERSH 0.318 0.043 7.443 0.000 0.318 0.496
## LEARNING 0.390 0.049 7.884 0.000 0.390 0.533
## ENVIRONM 0.368 0.049 7.571 0.000 0.368 0.506
## LEADERSH ~~
## LEARNING 0.386 0.047 8.259 0.000 0.386 0.566
## ENVIRONM 0.249 0.043 5.800 0.000 0.249 0.369
## LEARNING ~~
## ENVIRONM 0.332 0.050 6.626 0.000 0.332 0.430
## VALUE ~~
## TENURE 0.000 0.000 0.000
## JOBSTYLE ~~
## TENURE 0.000 0.000 0.000
## TEAMWORK ~~
## TENURE 0.000 0.000 0.000
## LEADERSH ~~
## TENURE 0.000 0.000 0.000
## LEARNING ~~
## TENURE 0.000 0.000 0.000
## ENVIRONM ~~
## TENURE 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .commit 1.456 0.123 11.853 0.000 1.456 0.586
## .out 0.291 0.025 11.853 0.000 0.291 0.739
## VALUE 0.647 0.055 11.853 0.000 0.647 1.000
## JOBSTYLE 0.907 0.077 11.853 0.000 0.907 1.000
## TEAMWORK 0.691 0.058 11.853 0.000 0.691 1.000
## LEADERSH 0.598 0.050 11.853 0.000 0.598 1.000
## LEARNING 0.776 0.065 11.853 0.000 0.776 1.000
## ENVIRONM 0.765 0.065 11.853 0.000 0.765 1.000
## TENURE 84.834 7.157 11.853 0.000 84.834 1.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## a1b 0.054 0.018 3.008 0.003 0.054 0.069
## a2b 0.030 0.014 2.186 0.029 0.030 0.046
## a3b 0.013 0.015 0.882 0.378 0.013 0.018
## a4b 0.014 0.015 0.898 0.369 0.014 0.017
## a5b 0.046 0.017 2.700 0.007 0.046 0.064
## a6b 0.035 0.015 2.333 0.020 0.035 0.048
## a7b 0.003 0.001 2.477 0.013 0.003 0.040
## total 0.014 0.004 3.869 0.000 0.014 0.206
学习成长到工作绩效的直接效应(d)为0.158,p值为0,说明该效应存在。比较前后模型的拟合度指标发现,卡方值明显下降,p值已大于0.05,RMSEA和SRMR都小于目标值。新模型完美拟合,无需再进行调整。
fitmeasures(path1_fit, c("chisq", "df", "pvalue", "nfi","nnfi","cfi", "rmsea", "srmr")) # 原设定模型
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 27.106 12.000 0.007 0.969 0.945 0.982 0.067 0.060
fitmeasures(path2_fit, c("chisq", "df", "pvalue", "nfi","nnfi","cfi", "rmsea", "srmr")) # 调整后模型
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 13.714 11.000 0.249 0.984 0.989 0.997 0.030 0.041
前后两个模型是嵌套关系,可以进行嵌套模型的卡方差异检验,检验模型修正带来的拟合优度增幅是否具有显著意义。
anova(path1_fit, path2_fit)
##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## path2_fit 11 7098.7 7222.4 13.714
## path1_fit 12 7110.1 7230.1 27.106 13.391 0.20999 1 0.0002528 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
可以绘制两个模型的路径图进行比较。
统合模型分析步骤:通常采用两阶段策略,第一个阶段是针对测量模型确定因素结构的拟合性,第二阶段则是在不改变测量模型的前提下,增加结构模型的设定,并评估结构模型的拟合性。该方法的优点在于测量模型的检测可以提供潜变量的聚敛与区分效度的信息,结构模型则可以提供预测效度的证据。
变量组合的策略:结构方程模型分析经常会以变量组合策略简化测量模型,使结构模型得以在比较简化的情形下进行估计。例如某个潜变量有6个测量变量,可以每两题加起来求平均值,将6个变量聚合成3个变量以降低模型的复杂程度。聚合层次的变量分数成为组合分数。特别是当样本量太少,而模型估计参数很多,变通做法是通过将测量变零组合成单一变量,将潜变量改为观察变量后进行路径分析。变量组合的优点包括:简化模型、提高模型拟合度、获得较理想的估计解、得到较好的测量信度、较高的变量解释力、观察变量的尺度更具有等距性与正态性、更理想的检验效力(每题的样本量与题数的比值更高)、避免特殊题项的干扰。但是,变量组合不应为简化而简化,需要提出简化的理由与统计量上支持,必须经过验证性因子分析的检验,保证构念的单一维度,并且组合的题项要具有相同的尺度。
250位教师的创意教学行为调查数据,包括教师的创造人格特质(Person)、创意教学自我效能感(Efficacy)、组织社会化(Socialized)、教学创新行为(Crea)等4个潜变量的10个测量变量。其中教学创新行为已将9个量表题项加总成组合分数。
研究假设: 1. 教师创意教学⾃我效能越⾼,越能够表现出创意教学⾏为。 1. 教师个⼈创造性格越强、组织社会化程度越⾼,则创意教学⾃我效能越⾼。 1. 教师个⼈创造性格与组织社会化程度两项特质,会透过⾃我效能的中介作⽤,间接影响教师的创意教学⾏为。
数据是以协方差矩阵的形式给出的。可以通过热力图展示测量变量间的相关系数。
library(tidyverse)
library(lavaan)
library(modelsummary)
library(semPlot)
library(corrplot)
## corrplot 0.95 loaded
library(semptools) # 给路径图标记显著性
N <- 250 # 样本量
# 协方差矩阵
COV <- structure(c(17.711, 1.843, 0.801, 1.243, 1.692, 7.518, 7.585,
6.499, 3.02, 2.663, 1.843, 0.404, 0.146, 0.227, 0.226, 1.074,
1.062, 0.906, 0.369, 0.318, 0.801, 0.146, 0.374, 0.11, 0.115,
0.301, 0.538, 0.388, 0.237, 0.307, 1.243, 0.227, 0.11, 0.589,
0.208, 0.531, 0.609, 0.434, 0.322, 0.335, 1.692, 0.226, 0.115,
0.208, 0.393, 0.741, 0.806, 0.664, 0.196, 0.326, 7.518, 1.074,
0.301, 0.531, 0.741, 7.565, 5.202, 4.53, 1.947, 2.092, 7.585,
1.062, 0.538, 0.609, 0.806, 5.202, 7.046, 4.626, 1.524, 1.824,
6.499, 0.906, 0.388, 0.434, 0.664, 4.53, 4.626, 6.335, 1.649,
1.662, 3.02, 0.369, 0.237, 0.322, 0.196, 1.947, 1.524, 1.649,
4.482, 2.335, 2.663, 0.318, 0.307, 0.335, 0.326, 2.092, 1.824,
1.662, 2.335, 2.982), dim = c(10L, 10L), dimnames = list(c("CREAT",
"SEFF1", "SEFF2", "SEFF3", "SEFF4", "PER1", "PER2", "PER3", "SOC1",
"SOC2"), c("CREAT", "SEFF1", "SEFF2", "SEFF3", "SEFF4", "PER1",
"PER2", "PER3", "SOC1", "SOC2")))
# X和Y相关系数,等于二者的协方差除以二者标准差的积:
COR <- COV / (sqrt(diag(COV) %*% t(diag(COV))))
# 调用corrplot命令,绘制热力图
# Generate a lighter palette
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(COR, method = "shade", shade.col = NA, tl.col = "black", tl.srt = 45,
col = col(200), addCoef.col = "black", cl.pos = "n")
模型拟合情况大致良好,卡方p值为0,卡方自由度比为2.6,接近2,NFI为0.94,CFI为0.96,RMSEA为0.08,SRMR为0.04,可以继续检验MI指数按照验证性因子分析介绍的方式调整测量模型,这里直接进行第二阶段的结构模型估计。
cfa_mod <- '
# measurement model
CREA =~ 1*CREAT # 由于CREA只有一个测量变量,因此默认测量变量残差方差为0
SEFF =~ SEFF1 + SEFF2 + SEFF3 + SEFF4
PER =~ PER1 + PER2 + PER3
SOC =~ SOC1 + SOC2'
cfa_fit <- cfa(model = cfa_mod, sample.cov = COV, sample.nobs = 250)
fitmeasures(cfa_fit)[c('chisq', 'df', 'pvalue', 'nfi', 'nnfi', 'cfi', 'rmsea', 'srmr')] |> round(3)
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 80.964 30.000 0.000 0.936 0.937 0.958 0.082 0.042
模型设定和拟合情况如下,可以看到,模型的拟合程度还是不错的,但仍有修正空间。
sem_mod<-'
# measurement model
CREA =~ 1*CREAT
SEFF =~ SEFF1 + SEFF2 + SEFF3 + SEFF4
PER =~ PER1 + PER2 + PER3
SOC =~ SOC1 + SOC2
# structural model
SEFF ~ a1*PER + c1*SOC
CREA ~ a2*PER + b1*SEFF + c2*SOC
#set the correlation between variables to free
SOC ~~ PER # 外生潜在变量间的协方差
#define indirect effect (a*b)
a1b1 := a1*b1
c1b1 := c1*b1'
sem_fit <- sem(model = sem_mod, sample.cov = COV, sample.nobs = 250)
fitmeasures(sem_fit)[c('chisq', 'df', 'pvalue', 'nfi', 'nnfi', 'cfi', 'rmsea', 'srmr')] |> round(3) # 查看拟合结果
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 80.964 30.000 0.000 0.936 0.937 0.958 0.082 0.042
p1 <- semPaths(sem_fit,
what = "std",
fade = F, # 关闭路径颜色渐变
rotation = 2,
layout = "tree2",
edge.color = "darkgrey",
esize = 3,
edge.label.cex = 1.1,
ask = FALSE)
p2 <- mark_sig(p1, sem_fit) # add asterisks
plot(p2)
模型修饰指数(MI)最高的五项。一般来说,基于MI修饰模型时,必须考虑到切实的理论意义。
从给出的建议看,可采纳一下几条(具体思考历程,请参见原书253-254页):
遵照第三行,增加CREA到SEFF4的测量因素载荷(原书253强调,不要轻易改动测量模型,因此,这里的改动只是出于示范目的)。 增加SEFF3和SEFF4的残差相关性。
modificationIndices(sem_fit) |>
select(1:5) |>
arrange(desc(mi)) |>
head(5) |>
knitr::kable(caption = "排名前五的模型修饰指数")
| lhs | op | rhs | mi | epc |
|---|---|---|---|---|
| PER | =~ | SEFF1 | 13.127181 | 0.1543965 |
| CREAT | ~~ | SEFF4 | 13.031960 | 0.3288959 |
| CREA | =~ | SEFF4 | 12.759405 | 0.0733377 |
| SEFF1 | ~~ | SEFF4 | 11.869383 | -0.0678181 |
| SEFF2 | ~~ | PER1 | 9.984666 | -0.1958852 |
拟合模型,并绘制路径图。
final_mod <- '
# measurement model
CREA =~ CREAT + SEFF4 # 修饰1
SEFF =~ SEFF1 + SEFF2 + SEFF3 + SEFF4
PER =~ PER1 + PER2 + PER3
SOC =~ SOC1 + SOC2
CREAT ~~ 0 * CREAT
# structural model
SEFF ~ a1*PER + c1*SOC
CREA ~ a2*PER + b1*SEFF + c2*SOC
#set the correlation between variables to free
SOC ~~ PER # 外源潜在变量间的协方差
SEFF3 ~~ SEFF4 # 修饰2
#define indirect effect (a*b)
PERindiCREA := a1*b1
SOCindiCREA := c1*b1
# define total effect
PERtotCREA := PERindiCREA + a2
SOCtotCREA := SOCindiCREA + c2
'
final_fit <- sem(model = final_mod, sample.cov = COV, sample.nobs = 250)
p1 <- semPaths(final_fit,
what = "std",
fade = F, # 关闭路径颜色渐变
rotation = 2,
layout = "tree2",
edge.color = "darkgrey",
esize = 3,
edge.label.cex = 1.1,
ask = FALSE)
p2 <- mark_sig(p1, final_fit) # add asterisks
plot(p2)
由于最终模型和原模型有嵌套关系,因此可以进行卡方检验,显示模型具有显著改进。
anova(final_fit, sem_fit)
##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## final_fit 28 7919.1 8014.2 59.089
## sem_fit 30 7937.0 8025.0 80.964 21.875 0.19937 2 1.778e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
查看最终模型的拟合结果:
fitmeasures(final_fit)[c('chisq', 'df', 'pvalue', 'nfi', 'nnfi', 'cfi', 'rmsea', 'srmr')] |> round(3) # 查看拟合结果
## chisq df pvalue nfi nnfi cfi rmsea srmr
## 59.089 28.000 0.001 0.953 0.959 0.974 0.067 0.039
展示创造性格PER和组织社会化SOC,对于创意教学行为CREA的间接效应和总效应量及其统计检验。
summary(final_fit, standard = T)$pe |>
filter(op == ":=") |>
select(label, est, se, z, pvalue, std.all) |>
mutate(across(where(is.numeric), ~round(., 3))) |>
arrange(label) |>
knitr::kable(caption = "潜在变量路径分析的部分复杂效应量")
| label | est | se | z | pvalue | std.all |
|---|---|---|---|---|---|
| PERindiCREA | 0.481 | 0.142 | 3.396 | 0.001 | 0.262 |
| PERtotCREA | 1.441 | 0.123 | 11.681 | 0.000 | 0.783 |
| SOCindiCREA | -0.053 | 0.064 | -0.835 | 0.404 | -0.019 |
| SOCtotCREA | -0.006 | 0.168 | -0.034 | 0.973 | -0.002 |
参考文献:《结构方程与建模的原理与应用》,邱皓政、林碧芳著,中国轻工业出版社,2009.